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_ . ABSTRACT 

. Context. The ionization of hydrogen in the solar chromosphere and transition region does not obey UTE or instantaneous statisti- 

' cal equilibrium because the timescale is long compared with important hydrodynamical timescales, especially of magneto-acoustic 

shocks. Since the pressure, temperature, and electron density depend sensitively on hydrogen ionization, numerical simulation of the 
solar atmosphere requires non-equilibrium treatment of all pertinent hydrogen transitions. The same holds for any diagnostic applica- 
tion employing hydrogen lines. 

Aims. To demonstrate the importance and to quantify the effects of non-equilibrium hydrogen ionization, both on the dynamical 
' (— I ' ' structure of the solar atmosphere and on hydrogen line formation, in particular Ho-. 

' Methods. We implement an algorithm to compute non-equilibrium hydrogen ionization and its coupling into the MHD equations 

I , within an existing radiation MHD code, and perform a two-dimensional simulation of the solar atmosphere from the convection zone 

Q ■ to the corona. 

V-^ ' Results. Analysis of the simulation results and comparison to a companion simulation assuming UTE shows that: a) Non-equilibrium 

I computation delivers much smaller variations of the chromospheric hydrogen ionization than for LTE. The ionization is smaller within 

, shocks but subsequently remains high in the cool intershock phases. As a result, the chromospheric temperature variations are much 

larger than for UTE because in non-equilibrium, hydrogen ionization is a less effective internal energy buffer. The actual shock tem- 
peratures are therefore higher and the intershock temperatures lower, b) The chromospheric populations of the hydrogen n = 2 level, 
which governs the opacity of Ka, are coupled to the ion populations. They are set by the high temperature in shocks and subsequently 
remain high in the cool intershock phases, c) The temperature structure and the hydrogen level populations differ much between the 
chromosphere above photospheric magnetic elements and above quiet internetwork, d) The hydrogen n = 2 population and column 
density are persistently high in dynamic fibrils, suggesting that these obtain their visibility from being optically thick in Ho- also at 
low temperature. 
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. 1. Introduction recombination timescales, the latter being far slower than the for- 
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The chromosphere repres ents the least understood regime of the | . j . | . | . 

sun (I.Tudge & Peteill998h . In this paper we address the treatment ICarlsson & Stein| iml 119951 120^ computed dynamical 



of hydrogen ionization in simulations of the solar chromosphere. ID simulations of the solar atmosphere including a detailed non- 
c3 ■ It is of paramount importance because hydrogen makes up 90% equilibrium treatment of the hydrogen rate equations including 
of the nuclei in the solar atmosphere, is an important and often ionization and recombination (i.e. not instantaneous statistical 
dominant electron donor, and contains a large part of the inter- equilibrium but employing pertinent time derivatives in the pop- 
nal energy of the gas. Hence, the ionization state of hydrogen elation rate equations, as in Appendix A of the present paper), 
strongly influences the temperature, pressure and electron den- I" the first paper they showed that the non-equilibrium eff^ects 
sity. Radiation magnetohydrodynamic (MHD) simulations of the lead to significant increase in shock temperature compared with 
atmosphere must therefore account properly for hydrogen ion- the case of instantaneous LTE ionization and recombination. In 
ization. This is not only important for the structure of the atmo- the second paper they supply a detailed analysis of the hydro- 
sphere within a simulation, but also for subsequent computation gen ionization in their simulations. They found that the chromo- 
of the emergent spectrum from the simulation for comparison spheric hydrogen ionization/recombination timescale is of the 
with observations order of 50 s within hot shocks and 10^ - 10^ s in the cool in- 
iKlein et al.l ([T9761 [T978h and ([T980I) already showed tershock regions, and that hydrogen becomes partially ionized 
from ideahzed one-dimensional (ID) models that the assump- ^i*in shocks but, owing to the long recombination timescale, 
tion of instantaneous statistical equilibrium (SE) for hydrogen ^oes not recombine in the subsequent post-shock phase. As a 
ionization does not hold in a dynamical atmosphere containing consequence, the degree of ionization of hydrogen in the chro- 
shock waves. The temperature diff^erence between the hot shocks mosphere is rather constant with time, in stark contrast to what 
and the cool intershock phases produces disparate ionization and *e classical assumptions of statistical equilibrium or LTE would 

predict. 
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The present limitations on computing power do not yet per- 
mit such full-fledged non-equilibrium treatment of hydrogen 
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ionization including radiative transfer in 2D and 3D simulation 
geometry. Therefore, approximations remain necessary to make 
the problem computation ally tractable. I n this work we employ 
the method formulated bv lSolluml(ll999l) . 

Solium showed that in ID dynamical modeling it is pos- 
sible to avoid detailed evaluation of the radiation field in the 
relevant hydrogen transitions by prescribing suitable radiative 
rates. He found that this approximation works well up to just 
below the transition region. Sollum's method is fast enough 
for use in multi-dimensional geometry. It was earlier imple- 
mented in the code CO^BOLD tFrevtag et al. 2002) to perform 
a 3D non- magnetic hydrodynamical simulation of the chromo- 
sphere by iLeenaarts & Wedemever-Bohml (l2006h. This simu- 
lation confir med the conclusion obtained by Carlsson & SteinI 
(II992ll2002h from ID simulation: the degree of hydrogen ion- 
ization in the middle and upper chromosphere is determined by 
the passage of high-temperature shocks, irrespective of the cool 
intershoc k phases. It is relatively constant at abou t lO^^* to 10"^. 
However. ILeenaarts & Wedemever-Bohml ( |2006|) did not imple- 
ment back-coupling of the non-equilibrium ionization into the 
equation of state for the gas within the hydrodynamical simu- 
lation. This is diflicult to do in the CO^BOLD code because it 
employs an approximate Riemann solver for the hydrodynam- 
ics. 

In this paper we first present our implementation of Sollum's 
method including back-coupling of the non-equilibrium ioniza- 
tio n into the equati on of state in the Oslo MHD code described 
bv lHansteenI (|2004|) which uses a finite difference scheme for the 
hydrodynamics. We then discuss an extensive 2D simulation of 
chromospheric fine-structure evolution with this code and ana- 
lyze the results in terms of the hydrogen ionization balance, with 
separation between the chromospheric behavior above a mag- 
netic element and in an area resembling quiet-sun internetwork. 
The effects of non-equilibrium hydrogen ionization are demon- 
strated quantitatively through comparison with a similar com- 
panion simulation. The latter was started with identical initial 
conditions but its hydrogen ionization was set to obey LTE, i.e. 
instantaneous Saha-Boltzmann partitioning. 



2. Method 

In this section we describe our implementa tion of Sollum' s 
method in the MHD code developed in Oslo (iHansteenl l2004t) . 
This code is based on the staggere d grid code descr ibed in 
iDorch & Nordlundl ([1998 ), Mackay & GalsgaardI ( 1200 ll) and by 
Galsgaard & Nor dluncQ. It includ es the multi-group opacity- 
binnin g method of iNordlundId 1 9821) and the treatment of scatter- 
ing of lskartlienl (|2000|) for radiative transfer in the photosphere 
and chromosphere. In the transition region and corona it em- 
ploys optically thin radiative cooling. In addition, it treats radia- 
tive losses in strong lines and continua of hydrogen and Ca II in 
the upper chromosphere and transition region using an approxi- 
m ation based on detailed ID computation with the RADYN code 
of lCarIsson& SteinI (II992h . Thermal conduction along magnetic 
field lines is taken into account {Hansteen 2004}). 

In this section we first summarize the temporal evolution 
scheme of the MHD algorithm and then specify the modifica- 
tions through which we insert non-equilibrium hydrogen ioniza- 
tion, the boundary conditions, and our simulation setup. 



2.1. MHD evolution scheme 

The MHD solver is t he explicit thir d-order predictor-corrector 
scheme developed by iHy mani d 1 9791) but modified to allow vari- 
able timesteps. Its fundamental variables are the density p, mo- 
mentum p, internal energy density ei and the magnetic field B. 
The evolution equations for the fundamental variables depend 
only on the variables themselves, the temperature, the gas pres- 
sure, and their spatial derivatives. These derivatives are com- 
puted using a sixth-order scheme. The temperature and pressure 
are looked up in precomputed tables as function of density and 
internal energy. This table is computed assuming LTE except 
where the density and internal energy have typical coronal val- 
ues. There coronal equilibrium is assumed. 

The evolution equation for fundamental variable / is 

where / only depends on quantities known at time f. Hyman's 
scheme solves this equation from timestep n to timestep n + 1 as 
follows: the predictor step is 

/,S (2) 

and the corrector step 

/«+ 1 = fl2/„- 1 + ( 1 - fl2)/« + b2f„ + C2/,|*\ , (3) 

where oi, 02, bi, b2 and C2 are coefficients that depend on the 
current and previous timestep sizes. 



2.2. Implementation of non-equilibrium hydrogen ionization 

In order to compute non-equilibrium hydrogen ionization one 
has to solve the hydrogen rate equations 

drii 
'dt 



-i+V-{,vv)^J]njPj^-niJ]Pij, (4) 

where n, is the population of hydrogen level ;, v the macroscopic 
velocity, n\ the number of levels and Pjj the transition rate co- 
efficient between levels ; and j. The left-hand side represents a 
continuity equation for the hydrogen populations, the right-hand 
side a source term describing the transitions between the hydro- 
gen levels. These equations are solved using operator splitting. 
The continuity part 



at 



(5) 



See |http : //www . astro . ku . dk/~k:g/| 



is solved using Hyman's scheme in tandem with the fundamen- 
tal variables. It is not possible to use sixth-order spatial deriva- 
tives for the hydrogen populations because negative populations 
then arise occasionally from the steep gradients in the population 
densities. Instead, a positive-definite first-order upwind scheme 
is used to ensure positivity of the populations. 

After the predictor step, the rate part of the equations, 

is integrated over the timestep At while enforcing charge conser- 
vation, hydrogen nuclei conservation, energy conservation and 
instantaneous chemical equilibrium between atomic and molec- 
ular hydrogen. This yields predicted values of the hydrogen pop- 
ulations, temperature and pressure. Subsequently, the corrector 
step is performed for the fundamental variables and the advec- 
tion part of the hydrogen populations. After this step the rate 
equations, charge, energy and particle conservation and chemi- 
cal equilibrium are solved again to obtain the hydrogen popula- 
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tions, temperature and pressure for the new timestep. The algo- 
rithm can be summarized as follows: 

- predict the fundamental variables and advection of hydrogen 
populations; 

- solve the rate equations and the conservation equations for 
the predicted temperature and pressure; 

- correct the fundamental variables and advection of hydrogen 
populations; 

- solve the rate equations and the conservation equations for 
the corrected hydrogen populations, temperature and pres- 
sure. 

The radiative rate coefficients Rij that enter in the total rate 
coefficients P,y are computed using Sollum's method. For each 
radiative transition a depth-dependent radiation temperature Trad 
is prescribed. It is set equal to the local gas temperature in 
the deep layers of the atmosphere, ensuring LTE populations 
there. Above a certain height the radiation temperature follows 
a smooth transition to a predefined value and then becomes con- 
stant at that value. The latter values were determined by requir- 
ing that Sollum's method reproduces the non-equilibrium hy- 
drogen populations obtained in ID modeling of the solar chro- 
mosphere with the RADYN code. All radiative Lyman transi- 
tions are set to obe y detailed balancing. Further detail is gi ven in 
ISollumI (Il999h and lLeenaarts & Wedemever-Bohml (l2006h . 

The pressure and temperature are explicitly computed in 
this modification; the equation of state tables are not used ex- 
cept at the lower and upper boundaries (Sec. 12.31 ). The radiative 
losses due to hydrogen are still computed time-independently 
using Skartlien's multi-group scheme which employs tabulated 
group-mean opacities, scattering probabilities and Planck func- 
tions based on LTE populations. 

The rate equations, the energy, charge and nucleus conser- 
vation equations and the chemical equilibrium equation together 
form a set of coupled non-linear equations that are solved using a 
Newton-Raphson scheme. These equations and their derivatives 
are specified in the Appendix. 

2.3. Boundary conditions 

The lower boundary condition enforces LTE hydrogen popula- 
tions by solving for Saha-Boltzmann equilibrium instead of the 
rate equations. The hydrodynamic conditions regulate the mass 
flow across the boundary. Consistency between the total hydro- 
gen number density and the mass density is automatically en- 
forced through the equation of hydrogen nucleus conservation. 
Thus, it is not necessary to specify hydrogen level population 
flows across the lower boundary, and the population continuity 
equations do not have to be solved in the grid cells at the lower 
boundary. 

The upper boundary, located in the corona, uses the same 
approach but here the non-equilibrium rate equations are solved 
instead of the Saha-Boltzmann equations. Consistency between 
the mass density and the hydrogen density is enforced by adding 
or removing protons. This boundary condition produces coronal 
equilibrium because the ionization relaxation timescale of 0. 1 s 
is small compared to the dynanu cal timescale in the corona (see 
Fig. 6 of lCarlsson & Steinll2002h . 

The equation of state tables are used for the boundary con- 
ditions on the hydrodynamic variables. The tabulated values are 
accurate at the boundaries, because they are based on LTE or 
coronal equiUbrium consistent with the local hydrogen popula- 
tions. 



2.4. Simulation setup 

We have performed a two-dimensional simulation including 
non-equilibrium hydrogen ionization in a setup similar to the 
one used by Hansteen et al. (2006]) and lDe Pontieu et al. (2007|). 
It has a horizontal extent of 16.64 Mm with a resolution of 
32.5 km and 512 cells in the horizontal direction. The vertical 
extent is II.I Mm with 150 cells. The vertical resolution varies 
from 32 km in the convection zone to 440 km in the corona. 
Continuum optical depth unity is located about 1 .5 Mm above 
the lower boundary. The horizontal boundary condition is peri- 
odical. Both the lower and upper boundary conditions are open, 
allowing flows to leave and enter the box. The upper bound- 
ary strives to maintain the temperature there at 800,000 K be- 
cause the current two-dimensional models do not have suffi- 
cient magnetic dissipati on to maintain a corona self-consistently 
(iDe Pontieu et al.ll2007h . We use a five-level plus continuum hy- 
drogen model atom to compute the non-equilibrium hydrogen 
ionization. A heating term is added to the energy equation that 
drives the temperature back to 2,400 K when it falls below that 
value. This value has no physical meaning but is present for sta- 
bility reasons. 

The simulation was started from a relaxed snapshot of a pre- 
vious simulation which employed LTE ionization and ran for 
30 min of solar time. The effects of the LTE ionization dis- 
appeared in approximately 5 niin of solar time (see Fig. 1 of 
iLeenaarts & Wedemever-Bohml (|2006|) ). 

3. Results 

Figure [T] shows a snapshot of the simulation. Panel a displays 
the temperature. Some magnetic field lines that extend into the 
corona have been overplotted. All coronal field lines are rooted 
in two photospheric magnetic field concentrations at x = 4 and 
,Y = 1 1 Mm. Henceforth we refer to these concentrations as mag- 
netic elements and to the areas between them as internetwork. 

The areas without field lines are not field free, as can be 
seen in panel b which shows the magnitude of the magnetic 
field strength. The temperature panel displays granulation at 
z = Mm, a shock-ridden chromosphere up to 2 - 4 Mm height, 
and a hot corona reaching peak temperatures of lO*" K. The 
height and shape of the transition region strongly depend on the 
magnetic field configuration, with the corona reaching deeper 
down above the magnetic elements. Panel c shows the mass den- 
sity. It reaches minimum values in the transition region above 
internetwork, which consists of extended high-rising arcs (black 
in the panel). 

Panel d shows the non-equilibrium hydrogen ionization de- 
gree. It has a minimum between and 0.5 Mm and smoothly 
increases upward to the completely ionized corona. The chro- 
mospheric ionization degree does not follow the local gas tem- 
perature. Panel e shows the fraction of hydrogen nuclei bound 
in H2 molecules, which peaks in cool chromospheric regions be- 
tween 0.5 and 2 Mm. The black-purple noise above z - 2 Mm 
is a numerical artifact caused by the single precision output of 
the code. The code itself uses double precision to avoid such ar- 
tifacts. Panel f shows the population density of hydrogen in the 
« = 2 level, the lower level of the Ho- line. It roughly follows 
the density structure, with the exception of the transition region 
where it shows a sharp increase at the locations of the arcs of 
minimum mass density. 

Panel g shows the non-LTE population departure coefficient 
of the ground state of hydrogen bi. The ground state population 
remains close to LTE throughout most of the photosphere and 
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Fig. 1. Snapshot cutouts from the simulation, showing various quantities in a vertical plane after 8.5 minutes of solar-time evolution. 
Panel a: gas temperature, with magnetic field lines that extend into the corona overplotted in black; b: magnetic field strength; 
c: mass density d: non-equilibrium ionization degree of hydrogen; e: fraction of hydrogen atoms in the form of H2 molecules; f: 
hydrogen n - 2 level population; g: departure coefficient for the hydrogen n - I level population; h: departure coefficient for the 
hydrogen n = 6 level population; i: departure coefficient for the hydrogen n-2 level population. The columns used in Figs.|2]and[3] 
are indicated by dotted lines in panel d. 



chromosphere, except in strong chromo spheric shocks where 
there is under-ionization compared to LTE. The ground state is 
strongly overpopulated in the transition region and is in coronal 
equilibrium in the corona. 

Panel h shows the departure coefficient of the hydrogen ion 
density b(,. It is much larger than unity in chromospheric cool 
intershock regions and smaller than unity within chromospheric 
shocks. This demonstrates that the non-equilibrium ionization 
degree is higher than in LTE in intershock areas and lower in 
shocks. 

The departure coefficient of the n-2 level in panel i shows 
the same structure in the photoshere and chromosphere as b(, due 
to the strong coupling of the continuum and the excited neutral 
levels. In the corona, b2 is around 5 x lO'*, its coronal equilibrium 
value. 

Figure |2] shows the behavior of atomic hydrogen along the 
two columns marked in Fig [1] These were selected to sample 
a magnetic element (left-hand column) and quiet internetwork 
(right-hand column). Panels a and b show the temperature and 
mass density. The corona starts much lower above the magnetic 
element than above the internetwork. Panel a shows a strong 
shock at 1 Mm which will become a dynamic fibril when it 



reaches the corona, pushing it upward. In contrast, internetwork 
panel b shows no strong shocks. The density at the transition 
region is much lower in this case. 

Panels c and d show the non-equilibrium degree of ion- 
ization of hydrogen as thick curves. It reaches a minimum of 
10"^ at around 0.5 Mm and increases smoothly towards com- 
plete ionization in the corona. The corresponding LTE ion- 
ization, obtained from the simulation temperature and elec- 
tron density stratifications with the Saha-Boltzmann equations, 
is shown as a thin curve. The dramatic differences between 
the curves demonstrate the failure of instantaneous LTE ion- 
ization in the chromosphere and transition region. In the non- 
equilibrium case, the slowness of ionization and recombina- 
tion prevents total ionization in the shocks and fuU recombina- 
tion in their wakes, producing far smoo ther ionization behavior 
with time than LTE would predict (s ee ICarlsson & Steinll2002l ; 
iLeenaarts & Wedemever-Bohml2006h . Note that the LTE curves 
reach complete ionization in the transition region at slightly 
lower heights than the non-equilibrium ionization curves. Since 
hydrogen becomes the dominant electron provider already at 
10 ionization, the electron density equals the proton density 
except in the temperature minimum. 
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Fig. 2. Properties of the simulation along a col- 
umn in a magnetic element (left-hand column) 
and in the internetwork (right-hand column). 
Panels a and b: temperature (solid) and mass 
density (dashed, right-hand scale); c and d: 
non-equilibrium (thick) and LTE ionization de- 
gree (thin); e and f: non-equilibrium (thick) and 
LTE proton density (thin); g and h: population 
of the n = 2 level for the non-equilibrium (thick) 
and LTE (thin) case. 



Panels e and f show the non-equilibrium (thick) and LTE 
(thin) proton densities. Since hydrogen becomes the dominant 
electron provider already at 10 ionization, the electron den- 
sity equals the proton density except in the ionization minimum 
around 0.5 Mm (panels c and d). 

Panels g and h show the population density of the n-2 level 
of hydrogen «2 (thick: non-equilibrium; thin: LTE). It shows the 
same trend as the proton population, because all excited hydro- 
gen levels are strongly coupled to the proton population reser- 
voir. Thus, the small temporal variation of the ionization is fol- 
lowed by the n-2 population, and, therefore, the Ha opacity. In 
particular, the peaks in no in the transition region in both panels, 
at z = 1.8 Mm and z - 4.2 Mm, respectively, provide significant 
Ha visibility. 

Figure [3] shows the temporal evolution of the temperature 
in a magnetic element and in the internetwork. The magnetic 
element acts as a wave guide, in which shocks travel upward 
with a period of three minutes. When they reach the corona they 
push it upward. With time, these motions produce characteristic 
parabolic height variation. The same behavior is observed in so- 
called Ho- dynamic fibrils. The rece nt papers bv lHansteen et al.l 
(|2006|) and lDePontieu et al.l (|2007|) show that Ho- dynamic fib- 
rils are the observational signature of such magneto-acoustic 
shocks. Notice that the material in the wake of the shocks can 
be as cool as 2,400 K. 

The internetwork temperature is less structured. A movie of 
the gas temperature (available in the online material for this ar- 
ticle) shows that the internetwork chromosphere is pervaded by 
shocks which originate in the photosphere but often travel side- 
ways, away from the magnetic elements, and sometimes even 
downward. The transition region is located higher than above 
a magnetic element and exhibits relatively small and irregular 
temporal variations in height. 



The second row shows the ionization degree of hydrogen. 
It is rather constant in time in the chromosphere. Above 1 Mm 
height, it stays at about 1 % ionization, both in the magnetic el- 
ement and the internetwork. The transition to coronal temper- 
atures is smoother in the internetwork and the increase in ion- 
ization is correspondingly smoother The third row displays the 
population of hydrogen in the « = 2 level «2- It is very low in the 
almost completely ionized corona. The transition region shows a 
local maximum, which is persistent in time. The « = 2 population 
is higher in the transition region of the magnetic element than in 
the internetwork. The width of the transition region maximum 
in the upward phase of the dynamic fibrils decreases suddenly 
when the fibril descends again. 

The fourth row displays the column density of «2, which is 
proportional to the vertical optical depth in the Ha line. In the 
internetwork, the column density starts to increase higher in the 
atmosphere than in the magnetic elements, because the chromo- 
sphere extends to larger heights in the internetwork. However, 
the column density at the top of the dynamic fibrils (the top 
of the light blue bulges) is lO'"* cm"^, two orders of magnitude 
larger than the internetwork column density at equal height. The 
same column density in the internetwork is reached at a height 
of 0.8 Mm. 

3.1. Comparison with companion LTE simulation 

Figure |4] presents a comparison between the simulation with 
non-equilibrium hydrogen ionization and a simulation with the 
old code which employs LTE ionization. Both simulations were 
started at r = min from the same relaxed snapshot computed 
with LTE ionization. 

The bottom panels of|4]show that at height z= 1 Mm the dif- 
ferent treatments of ionization produce only a slight difference 
in temperature variation. The wave pattern is almost identical. 
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Fig. 3. Time slices of the gas temperature (first row), the ionization degree of hydrogen (second row), fraction of hydrogen in the n-2 
level (third row), and the n = 2 column density (fourth row) in a magnetic element (left-hand column) and in the internetwork (right- 
hand column). The upper-left magnetic element panel shows dynamic fibrils pushing the corona upward with 3 min periodicity. The 
upper-right internetwork panel shows rather unstructured shocks and a slowly varying height of the transition region. The snapshot 
used in Fig.[T]and|2]is indicated by a black dotted hne. 



but with non-equilibrium ionization the fluctuations are larger. 
At larger heights, the differences between the simulations are 
large (upper two rows). The shock patterns (thin hot threads) 
differ markedly. The shock temperatures are much larger in the 
non-equilibrium simulation. They typically are 10,000 K and 
sometimes reach 13,000 K, whereas the LTE shock temperatures 
do not exceed 8,000 K. These high temperatures are due to the 
comparative lack of ionization. Because the internal energy of 
the gas is not stored as hydrogen ionization energy, it remains 
as kinetic energy of the gas particles, raising the temperature 
(ICarlsson & Steinlll992h . 



The temperatures also differ greatly in the cool intershock 
regions. In LTE, they are between 3,000 and 5,000 K whereas 
the non-equilibrium simulation reaches intershock temperatures 
of about 2,500 K. These low values result from the reverse pro- 
cess: over-ionization compared with LTE. More energy remains 
stored as ionization energy, leaving less kinetic energy for the 
gas particles. 
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Fig. 4. Time slices of the gas temperature at 2, 1.5 and 1 Mm height (top to bottom). Left-hand column; simulation with LTE 
ionization; right-hand column: with non-equilibrium ionization. Both simulations were started with the same snapshot at r = min. 
The intermittently present red patches are caused by the corona, whose lower boundary is moved up and down by shocks traveling 
upward above the magnetic elements. The snapshot used in Fig.[T]is indicated with the horizontal dotted line; the columns used in 
Figs.|2]and[3]with vertical dotted lines. 

lation. It is not clear how low the actual chromospheric minima 
may reach because radiative heating in the hydrogen continua 
and other strong spectral features is not taken into account in the 
simulation, only their radiative cooling. 

4.2. Discussion 



From the analysis of our simulation we obtain the following 
picture. The internetwork chromosphere is irregularly pervaded 
by shocks. The temperature is typically 10,000 K in shocks 
and can be lower than 2,500 K in cool intershock areas. The 
transition region is arc-shaped, with the arc footpoints situated 
above the magnetic elements, and reaches a maximum height of 
4 Mm. The chromosphere above magnetic elements shows up- 
ward propagating shocks at about 3 minute periodicity, which 
push the corona upward from 1 .5 Mm to 3 Mm height. 

The chromospheric hydrogen ionization degree increases 
smoothly from 10"^ at 0.5 Mm to complete ionization in the 
corona and does not respond to temperature changes because of 
the slow recombination behind shocks. The hydrogen « = 2 pop- 
ulation in the chromosphere is coupled to the proton population 
and, as a consequence, varies only weakly with time too. The 



4. Discussion & Conclusions 

4.1. Limitations of the simulation 

Our implementation of non-equilibrium hydrogen ionization has 
various limitations. 

First, the assumption that all Lyman transitions ar e in de- 
tailed balance is justified up to the transition region dSolluml 
119991) . However, the transition region is optically thin in most 
Lyman features, requiring detailed radiative transfer modeling 
to evaluate their influence on the hydrogen populations. 

Second, the multi-group radiative transfer within the simula- 
tion, which sets the radiative cooling and heating, employs LTE 
ionization. For given internal energy and mass density, the radia- 
tive transfer uses the group-mean opacity, scattering probability 
and Planck function based on the corresponding LTE (or coro- 
nal equilibrium) temperature and electron density. The radiative 
cooling in the chromosphere and transition region, where devi- 
ations from equilibrium are largest, is thus inconsistent with the 
non-equilibrium temperature and electron density as computed 
in the simulation. 

Third, the cool parts of the simulation chromosphere often 
reach the limiting temperature of 2,400 K allowed in the simu- 
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n-2 population in dynamic fibrils above the magnetic elements 
is higher than in the internetwork chromosphere. The n-2 col- 
umn density, a measure of the optical depth of the Ha line, is 
two orders of magnitude larger in dynamic fibrils than at equal 
height in the internetwork. The same column density as at the top 
of the fibrils is only reached at a smaller height, about 0.8 Mm, 
in the internetwork. 

In the observatio ns of iHansteen et al. I (I2006h and 
iDe Pontieu et alj ( l2007h dynamic fibrils appear as dark "fingers" 
extending and retracting on top of a brighter background. 
They seem to be optically thick. Combining this observational 
appearance with the results of our simulation suggests that 
dynamic fibrils are optically thick in the Ha line core, whereas 
optical depth unity is reached far deeper in the internetwork 
atmosphere (equal column density in the lower panels of Fig.O 
so that the internetwork chromosphere adjacent to dynamic 
fibrils is optically thin. This explains that dynamic fibrils are 
observable along a slanted line of sight through the internetwork 
chromosphere. We also suggest that their low temperature com- 
bined with the effects of NLTE resonance scattering produces 
their low emergent intensity. The bright background is formed 
much deeper in the internetwork, where the higher temperature 
leads to higher emergent intensity. In summary, slow recombi- 
nation and strong coupling of the n-2 population to the ion 
population makes these fibrils opaque and therefore visible 
even when they are cool, regardless of the large n—2 excitation 
energy. Ha interpretation assuming a static atmosphere with 
instantaneous ionization and recombination, as frequently done 
in cloud modeling and inversions based on cloud modeling, is 
likely to be erroneous for such dynamic structures. Obviously, 
detailed radiative transfer computation based on time-dependent 
simulation with non-equilibrium ionization as done here is 
needed to properly assess Ha formation in dynamic fibrils. 

4.3. Conclusions 

We have presented an algorithm to compute non-equilibrium hy- 
drogen ionization with back-coupling to the equation of state in 
multidimensional radiation MHD simulations of the solar atmo- 
sphere. We performed a 2D simulation from the convection zone 
to the corona that employed this algorithm. From its analysis and 
comparison with a companion LTE simulation we conclude the 
following: 

- Inclusion of non-equilibrium hydrogen ionization is essential 
in simulations of the solar atmosphere because the resulting 
temperature structure and hydrogen populations differ dra- 
matically from their LTE values. 

- The degree of ionization of hydrogen in the chromosphere 
does not follow the local tem perature, as described already by 
ICarlsson & SteinI ( |2002[) and lLeenaarts & Wedemever-Bohml 
(120061) . Hydrogen is partially ionized in shocks but does not 
recombine in the cool shock wakes, owing to the slow recom- 
bination rate at low temperature (Figs. [11121 and|3]l. 

- Non-equilibrium hydrogen ionization causes more profound 
temperature variations in the chromosphere than would occur 
if LTE were valid (Fig.lU. The shock temperatures are higher 
and the intershock temperatures are lower, caused by the in- 
sensitivity of the hydrogen ionization degree to variations of 
the state parameters of the gas. 

- The population of the hydrogen n-2 level in the chromo- 
sphere is strongly coupled to the ion population. It therefore 
behaves approximately as the latter Its value is set by the 
high shock temperatures and subsequently remains high in 



the cool aftershock phases (Fig.|2|. This is quite contrary to 
the LTE prediction. 

- The simulation shows large differences in temperature struc- 
ture and hydrogen level populations between magnetic ele- 
ments and internetwork (Figs.|2land|3]l. 

- The Ha line opacity is proportional to the n-2 level pop- 
ulation; the Ha optical depth scales with the n-2 column 
density. Both are appreciably larger in dynamic fibrils than in 
the internetwork chromosphere at equal height (Fig. |3]l. We 
suggest that dynamic fibrils are optically thick in Ha and that 
their low temperature combined with scattering make them 
appear dark against the deeper-formed bright internetwork 
background. 

The next step is to compute Ha in detail from this simulation. 
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Appendix A: The equations and their derivatives 

The solution scheme requires partial derivatives with respect 
to the independe nt variables. The f ormul ation is similar to the 
RADYN code of ICarlsson feSteinl (IT9921) but since no detailed 
specification has been published so far, we supply the complete 
derivative list here. 

The present code solves the equations of chemical equi- 
librium, charge, energy and hydrogen nucleus conservation to- 
gether with the rate equations or with the Saha-Boltzmann equi- 
librium equations. These equations depend on the gas tempera- 
ture r, the electron density ng, the hydrogen population densi- 
ties «,, where «i is the ground state, «2-5 the four excited states 
and «6 the ionized hydrogen (proton) density, and the molecu- 
lar hydrogen density «h2- These are ten equations and nine un- 
knowns. We use only five of the six rate equations and eliminate 
the molecular hydrogen density, reducing the number of equa- 
tions and variables to eight. The equations are labeled Fi with 
/ = 1,...,8. 

The physical constants have their usual meaning: h, k, c, e 
and me are Planck's constant, Boltzmann's constant, the velocity 
of light, the electron charge and the electron mass. 
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As the mass density of the gas is known, the total number of 
hydrogen nuclei, whether in the form of H, H2 or bare protons is 
known and given by 

«g'=p«Hpei-g (A.l) 

where «Hperg is the number of hydrogen nuclei per gram solar 
gas. We denote the number density of nuclei of elements other 
than hydrogen as /inoH = "nuclei per h where «„ucieipeiH is the 
number of nuclei from elements other than hydrogen per hy- 
drogen nucleus. The number of electrons due to elements other 
than hydrogen per hydrogen nucleus is n"°". The internal en- 
ergy of elements other than hydrogen per hydrogen nucleus is 
gnoH- These quantities and their derivatives are determined nu- 
merically from tables as function of T and «e- The molecular hy- 
drogen density is nH2 and the internal energy of H2 per molecule 
is eH2- The latter is computed using polynomial approximations 
bv I Vardval (119651) 



Chemical equilibrium. The equation of chemical equilibrium 
between neutral atomic hydrogen and molecular hydrogen is 
given by 



K 



(A.2) 



where K is the chemical equilibrium constant whose value and 
its temperature derivativ e are given b y polynomial approxima- 
tions as function of T by lTsujil (Il973h . The molecular hydrogen 
density depends on T and «, {i- 1, . . . , 5). The derivatives are 



dn 



H2 



dT 

dnm 
dtii 



2(i:ii»/) 

K 



^' dK 
df' 



Charge conservation. The electron density is given by 



(A.3) 



(A.4) 



(A.5) 



which we rewrite in a form suitable for Newton-Raphson itera- 
tion: 



(A.6) 



tives are 






dFi 






drie 


111 




dFi 


1 




drit, 


We' 




dFi 








lie dT 





(A.7) 
(A.8) 
(A.9) 



with Xi the energy of hydrogen level /. This functional depends 
on ris, tii and T. The partial derivatives are 



dF2 _ I (3kT 
Oris ei\ 2 

dF2 _ _ 1 l3kT 
dill e\ \ 2 



,de. 



noH 



H dn 
dnm 



dni 



+ 1 



dnm 
+ em-^ — +X> 
dni 



(A. 11) 
(A. 12) 



dFo 
'dT 



3k 

y 



dnm 

"e + "noH + nm + T^-r + 
dT 



tott^^^noH dem^ dnm 
+n'S^:zr- + rim^^ + em 



dT 



dT 



dT 



(A. 13) 



Particle conservation. The conservation of the total number of 
hydrogen nuclei «™ is given by: 



= 1 - 4t 

"h 



( 6 



Z rii - 2ny 



= 0. 



(A. 14) 



n 



dni 



This functional depends on «, and T . The partial derivatives are: 

dFi 
dni 
dF3 

dn^ n 
dF3 
'dT 



l,(,.2f^|,,'=l,...,5, 

1 

tot' 
H 

2 dnm 



n'°' dT 



(A. 15) 
(A. 16) 
(A. 17) 



LTE populations. LTE is imposed to start up a new computation 
and at the lower boundary. The LTE hydrogen populations then 
obey the Saha-Boltzmann equations (/ = 1, . . . , 5) 

^3,, . 1 _ 1^2g,/2.»ar\^^^^_^^_^,,,, ^ ^^^^^ 

"e "6 gi \ I 

where gi is the statistical weight of level /. These equations de- 
pend on ng, n,, and T. Defining 



«e «6 gi \ j 

the partial derivatives are given by 



The functional Fi depends on n^, ng and T. The partial deriva- 



3+/ 



dnf, 
dFj^i 

dni 
dFi+i 

dne 
dF^ 

dT 



rii ' 

ne ' 
= -Ki 



Xb -Xi 
kT^ 



3 

TT 



(A. 19) 

(A.20) 
(A.21) 
(A.22) 
(A.23) 



Internal energy. The second functional specifies the distribu- Rate equations. Outside LTE, the change of the hydrogen pop- 
tion of the internal energy over the various contributions: 



Fi 



1-1 



3kT 



tie + HnoH + "02 + ^ «; 

1=1 

6 \ 



+ «H ^noH + nm em + 



Y^niXi 



= 0, 



(A. 10) 



ulations n\ . . . ng over a timestep Af can be expressed as 



riiitQ + At) - mito) = I Yj"j{^Ji + '^ji) 

6 

-rii {^U + Cij) df, 



(A.24) 
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with Rij and C,j the radiative and collisional rate coefficients 
from level / to level j. Defining 

(A.25) 
(A.26) 



Rij + Cij, 



and writing «, = n, (fo+Af) and n° - n, (fo) yields the discretized, 
implicit form of Eq. IA.24I 



Af 



( 6 



n n 



-1=0. 



(A.27) 



The equations are dependent on rij (j + i), rig and T. The 
partial derivatives are 



dF 



3+i 



1 -AtPii 



drii 
drij 

drie 
dT 



n1 
AtP 



ji 



n 
At 

Af 

Jo 



to 



dpj.\ 



(A.28) 
(A.29) 

(A.30) 
(A.31) 



Only five of the six equations I A. 241 are used to avoid overspeci- 
fication through ignoring the one with the largest value of n'". 

LIE population ratios. Evaluation of the rate coefficients in- 
volving the continuum level «6 often involve the LTE population 
ratio (/ = 1 , . . . , 5) given by 



gi I InmJiT 



-3/2 



with the derivatives 
dT \n* 



kT^ ^ 2T j nl ' 



(A.32) 

(A.33) 
(A.34) 



Collisional rate coefficients. The expressions for the collisional 
rate coefficients contain temperature-dependent coefficients Cexc 
and Cion for excitation and ionization. Th eir values and d eriva- 
tives are determined from a table based on lJohnsonI ( Il972h . 

The collisional rate coefficients for bound-bound transitions 
are, with j > i and Cij the upward rate coefficient: 



Cji = — HeCexcVr, 
8J 



i^ij — "e^exc V ^ 



with derivatives 
dCji^ ^ Cji_ 
drie rie 

^^Ji 8i ,r^l "^Cexc Cexc 

- —n^yT\ + 



dT 8J 
9Cij ^ Cij 
drie rie ' 



dT 2T 



(A.35) 
(A.36) 

(A.37) 
(A.38) 
(A.39) 



dCij 
dT 



\ dT 



+ Ce: 



Xj -Xi 



kT^ 



1 

TT 



I (A.40) 



For bound-free transitions the collision rate coefficients are, 
with C,6 the upward rate coefficient: 

C/6 = «eCion(r) Vf e-^'-^')/*^ (A.41) 
Cfi/ = lc,-6, (A.42) 



with derivatives 



dCi(, 
drie 

dCi6 
dT 

dCg 
5«e 
dCg 

dT 



Ci(, 

«eVre-^^'-^')/*^(^ + Ci, 



= 2^, 

dCi6 

dT 



if 



dT 



X6 -Xi 



Xb -Xi 



kT^ 



1 

TT 



kT^ 



(A.43) 
I (A.44) 
(A.45) 
(A.46) 



Radiative rate coefficients. The derivation of the radiative rate 
coefficients can be found in ISollum(1999) . The rate coefficients 
for bound-bound transitions are, with j > i and Rij the upward 
rate coefficient: 

4nV 2hvl 1 

(A.47) 



^'7 - -, /lu 



QhvalkT,-.^oi _ \ ' 



R 



2„2 



-/lu- 



2hvi 



1 



(A.48) 



■'' gj hvotneC''"' 1 - Q-hyolkT,.^i ' 

where /lu, vq and Trad are the oscillator strength, the line cen- 
ter frequency and the prescribed radiation temperature. In areas 
where Tjad = T, i.e., deep in the atmosphere, the temperature 
derivatives are: 



dRij _ 2,n^e~hvl 



dRi 



gi SnVhv^^ 



fij 



Jivo/kT 



(A.49) 



(A.50) 



dT gj m^c^kT^-"\e'"'o"^T _ 1)2' 

In areas where Tjad T, but constant, the temperature deriva- 
tives are zero. 

The upward radiative rate coefficient for bound-free transi- 
tions is: 



8;r 



^/6 = ^Q'oVo2^£il«^y, 

«=1 ' 



hvo 



rad 



(A.51) 



where q-q is the radiative absorption cross section at the ioniza- 
tion edge frequency vo and Ei{x) the first exponential integral 
with argument x. Its temperature derivative is: 



dRi6 Stt 3 1 1 



T — - -^onvr. , where Trad = r, 

dT T e'''^o/kT _ I ma , 

dRi6 



dT 



- where Tj-ad T. 
The downward bound-free radiative rate coefficient is: 



87r , " • v-i 
R6i ^ -^aovo-^ 2^Ei 

6 n=\ 

The derivatives are: 
dRei _ R6i 

dR6i ^ dRi6 «; 
dT dT nl 
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n h 1 

Tr-dd 



livo 
kT 



(A.52) 

(A.53) 
3: 

(A.54) 
(A.55) 



^6/1^ + ^1 where r,ad = T, (A.56) 
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— Q'ovo 



^ nl 



iT + T, 



(3 hvo \ 
^ + ^ I where T,^d + T. 



(A.57) 



